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ABSTRACT 

Gamma-ray bursts (GRBs) by virtue of their high luminosities can be detected up to very 
high redshifts and therefore can be excellent probes of the early universe. This task is hampered 
by the fact that most of their characteristics have a broad range so that we first need to obtain an 
accurate description of the distribution of these characteristics, and specially, their cosmological 
evolution. We use a sample of about 200 Swift long GRBs with known redshift to determine 
the luminosity and formation rate evolutions and the general shape of the luminosity function. 

In contrast to most other forward fitting methods of treating this problem we use the Efron 
Petrosian methods which allow a non-parametric determination of above quantities. We find a 
relatively strong luminosity evolution, a luminosity function that can be fitted to a broken power 
law, and an unusually high rate of formation rate at low redshifts, a rate more than one order 
of magnitude higher than the star formation rate (SFR). On the other hand, our results seem to 
agree with the almost constant SFR in redshifts 1 to 3 and the decline above this redshift. 

Subject headings: Gamma rays: bursts-cosmology: early universe-stars: formation-general meth¬ 
ods: statistical 


1. Introduction 

Observations of increasing numbers of gamma-ray bursts (GRBs) with measured redshifts, up to z ~ 10, 
by instruments on board BeppoSAX, HETE, INTEGRAL , and in particular Swift, have stimulated many uses 
of them as cosmological probes either as “standard Candles” (SC) for determination of global cosmological 
parameters, such as density parameters and equations of state, or as probes of the early phases of the universe 
such as reionization, or star formation rate (SFR) and cosmic metalicity evolution (CME) at high redshifts. 
Unfortunately, most intrinsic distance dependent charcteristics of GRBs, such as their peak luminosity, total 
emitted energy, etc, have a broad distribution making them unsuitable as SC. However, there has been 
several attempts to discover some distance dependent characteristic that shows a well defined correlation 
with another distance independent characteristic, which can then be used to determine distances as in Cepheid 
variables or type la supernovae. Example of such relations for GRBs are the correlations between lag and 
luminosity (Norris et al. 2000), the variability and luminosity (Reichart et al. 2001), and the peak energy 
E p of the vF v spectrum and the total (isotropic) gamma-ray energy £i so (Amati et al. 2002; Lamb et al. 
2004; Attiea et al. 2004) or the beaming corrected energy (Ghirlanda et al. 2004). For a more general 
review see Xiao & Schaffer (2009). 

There are, however, many uncertainties in the claimed relations rendering the cosmological tests unre¬ 
liable. First, it is now clear that these characteristics have broad distributions and that the correlations are 
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not simple one-to-one relations that is sometimes claimed and do not have the small dispersion required for 
precision cosmological tests. Instead, most of the correlations are statistical in nature, as originally predicted 
by Lloyd et al. (2000) long before any redshifts were measured, and not valid for the GRB population as a 
whole (Nakar & Piran 2004; Band & Preece 2005). Butler et al. (2009) stress that a careful accounting of 
observational selection effects is required, and suggest using the methods developed by Efron and Petrosian 
(EP for short, see below), to quantify the nature and degree of the correlations before any cosmological tests 
can be affected. 


Secondly, even if there exists a one-to-one relation; e.g. £i so = £qC(E p / Eq), the relation between the 
cosmological parameters, such as matter and dark energy density parameters f l m and Qa and observed 
quantities, such as redshift z, flux f(t) [or fluence F = J f(t)dt\ and E p , are complicated. Here £i so is related 
to the total gamma-ray energy fluence Ftot as 


£ iso = ^d 2 L F to t/Z, where d L 


{c/H 0 )Z [ dz'/y/nfi) 
Jo 


(1) 


is the luminosity distance, and f i(z) = e(z)/e o describes the evolution of the total energy density e(z), of all 
substance. 1 In this case the cosmological parameters are related to the observations as 


which involves at least 2 unknown functions f l(z) and C(E p ), and 3 more if the correlation function C, and 
the scales £q and Eq evolve in time. Without a knowledge of the forms of the evolution functions C[E p ,z ), 
£o(z) and E 0 (z) precise determination of the H(z) is not possible. Most attempts to this end have assumed 
not only that the correlation function has small dispersion (e.g. C(E p ) oc 5(E p — Eq)) but also that it remains 
narrow at all redshifts (i.e. Eq does not evolve), and that there is no luminosity or £ ISO evolution (i.e. £o is 
a constant). 

There has been some attempts to determine some aspects of the evolution (see e.g Li 2007), but most 
of these evolutionary trends have not been addressed. These evolutions, in principle, can be determined for 
an assumed cosmological model given a large enough sample. However, using such results to determine the 
cosmological model will be a circular and meaningless exercise. Thus, at this state of our knowledge, the use 
of the current GRB data for determining the global cosmological parameters seems premature. We need to 
learn more about the nature of the GRBs and the cosmological evolution of their characteristics before they 
can be used for this task. The immediate situation may be more akin to star forming galaxies and active 
galactic nuclei (AGNs), whose characteristics have also broad distributions and may evolve in time. This 
fact has shifted the focus of activity (both in galaxies, AGNs and GRBs) to the investigation of structure 
formation, the building process of the black holes, and to SFR and CME. 


In recent years there has been increased activity in trying to use the existing GRB data to determine the 
shape and evolution of the luminosity function (LF), 4>(L,z), the GRB formation rate density (co-moving) 
p(z) (GRBFR for short), and its relation to the SFR or CME (see e.g. Porciani & Madau 2001; Natarajan 
et al. 2005; Daigne et al. 2006; Jakobson et al. 2006; Le & Dermer 2007; Salvaterra et al. 2009; Butler 
et al. 2010; Wanderman & Piran 2010; Howell et al. 2014). All these works use the forward fitting (FF) 
method, whereby the observed data, such as flux, fluence and redshift distributions are fit to prediction 


1 Here Z = 1 + z, eo = 3 -Hqc 2 /(8ttG) and il(z) = £l m Z 3 + Pa- In what follows we use a Hubble constant Ho = 70 km/(s 

Mpc), and assume a flat universe with Q m = 0.3 for matter and Q,\ = 0.7 for the cosmological constant. 
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of models with some assumed parametric forms for the numerous functions [LF, p(z ), spectrum] and their 
evolutions. For example, a common practice is to ignore luminosity evolution (see, however, Salvaterra et 
al. 2012) and assume a GRBFR similar to the SFR, a power-law LF with breaks, Band spectrum with 
unique values of the high and low energy indexes, etc. A more objective procedure would be to determine 
these characteristics from the data directly and, as much as possible, non-parametrically. The statistical 
methods developed by Efron and Petrosian are designed exactly for this kind of analysis and have been used 
for analysis of cosmological evolutions of quasars (Maloney & Petrosian, 1999; Singal et al. 2011; Singal et 
al. 2013), GRBs (Lloyd et al. 1999); Kocevski & Liang 2006; Dainotti et al. 2013) and blazars (Singal et al. 
2012; Singal et al. 2014). 2 

In this paper we apply these method to data from Swift on long GRB with known redshifts with the 
aim of determining the cosmological evolution of the general LF which means the luminosity and formation 
rate variation with redshift. This parallels very closely with an earlier work using pre and early Swift data 
which formed a chapter of Aurelien Bouvier’s PhD thesis at Stanford (Bouvier 2010). The results of this 
work can also be found in Petrosian et al. (2013). In the next section we describe the data we use and in §3 
we present a brief description of the method as applied to these data. The results are presented in §4 and a 
brief summary and discussion is given in §5. 


2. Swift Data and Selection Effects 


Over the last decade, the Burst Alert Telescope (BAT) on board the Swift satellite has detected more 
than 800 long GRBs with peak fluxes in the 15-150 keV energy band and above the threshold flux / Pi ii m ~ 
2 x 10 -8 erg/(s cm 2 ). Approximately, 90% by the X-ray Telescope (XRT) and about one third by both 
XRT and the Ultraviolet/Optical Telescope (UVOT). Redshifts obtained for over 250 Swift GRBs from these 
instruments and the follow up observations they enable on larger, ground-based telescopes. Fig. |T] shows 
scatter diagram of luminosities and redshifts for 253 long GRBs, taken from Nysewander et al. (2007), 
Butler et al. (2009), and NASA’s online burst archive, which collates data from Evans et al. (2009) and the 
Gamma-ray Coordination Network circulars. 

From the redshifts and the observed peak gamma-ray energy flux / p , integrated over the observed Swift- 
band (15-150 keV), we calculate the peak luminosity for the assumed cosmological model and K-corrected 
for the same rest frame energy band as 3 


L = Attc1 2 l (z, D) f/K with K(Z) 


l™eV V E f(E) dE ' 


( 3 ) 


The spectra f(E) are fitted to either a power law f(E) = A(y^)“, a power law with an exponential cutoff 


2 Codes for application of these methods 
cran.r-proj ect.org/web/packages/DTDA/DTDA.pdf 


can be found at http://www.inside-r.org/node/99623 and 


3 Note that this definition of the K-correction follows the original definition that can be found in e.g. Peebles (1993). 
Sometimes the inverse of this is use relating a fixed rest frame energy band to a variable observed band, e.g Bloom et al. (2001) 

rl50 keV ^^(£/) d.E 

in which case K(Z) = 150 klv/z -• For a power law spectrum these give identical results but could be different when 

il5keV/Z Ef(E)dE 

there are significant spectral deviation from a simple power law. 
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(2+c)E 

f{E) = A{4 s ) a e-^~, 


or the Band model 


lf(E) 


( 2 +c)E 

A(j M)“ e Ep E K Ebr 

E > E br 


(a-P)E p 

2+Q' 


( 4 ) 


In Figure [T] the dotted and solid (black) curves show the truncation L > L min (z) = And 2 L (z,fl)f ln i n /K 
due to two fiducial peak flux thresholds flim and for K correction calculated for an average spectrum. To 
insure a higher completion level in what follows we use he larger limiting flux of /n m = 2 x 1CU 8 erg/(s cm 2 ) 
that includes 207 sources. The smoothed redshift distribution of this subsample is shown by the dotted (red) 
curve in the left panel of Fig. [4] The truncation induces a strong luminosity-redshift correlation, L oc z 3 - 3 , 
shown by the dashed (green) line, which complicates the determination of the true intrinsic correlation, 
namely the luminosity evolution. 

However, in addition to the bias involved in the detection of the prompt emission, the measurement of 
the redshift can introduce further truncations. X-ray and optical afterglows are vital for timely and accurate 
burst localization and redshift measurement. It is difficult to quantify the optical selection criterion. But as 
discussed in §4.3, we can include in our analysis the effects of X-ray selections. Since the optical and X-ray 
afterglow fluxes show relatively strong correlation, inclusion of X-ray selection effects may to some degree 
alleviate the problem of the optical selection effects. 


3. Methods and Approach 

The problem under consideration here requires determination of multi variate distribution from data 
truncated by observational selection effects. We will demonstrate our approach for determination of single 
LF and its evolution; z) from a flux limited sample such as that shown in Fig. [Q Without loss of 
generality, we can write the LF as 

L,z ) = p{z)ijj[L/g(z),a i ]/g{z), (5) 

where, in addition to the GRBFR p(z), we have introduced g(z) and a, to describe luminosity and shape 
evolutions, respectively. Given sufficiently large sample one can determine all three evolutions. Our experi¬ 
ence with evolution of the LF in other extragalactic sources such as quasars and blazars indicates that the 
least variable among these is the shape parameter(s). Thus, because of the limitations of the current GRB 
data we will focus here only on the GRBFR and luminosity evolutions and assume constant shape param¬ 
eters. The truncation or bias introduced by the flux limit [/ > /n m or L > L min (z) = ATrflj^(z)fn m /K(z)] 
is known as the Malmquist Bias. Many papers have been written to correct for this bias in astronomical 
literature (Malmquist 1925); Eddington 1940; Trumpler & Weaver 1953; Neymann & Scott 1959). Most of 
these early procedures assumed simple parametric forms (e.g Gaussian) distributions. However, since the 
discovery of quasars the tendency has been to use non-parametric methods: e.g (U/Umax) (Schmidt 1972; 
Petrosian 1973) or the C~ (Lynden-Bell 1971). For a more detailed description see review by Petrosian 
(1992). All these methods however, suffer from a major shortcoming because they assume that the variables 
are independent or uncorrelated, or that the LF is separable; ^(L, z) = ip(L)p(z). This ignores luminosity 
evolution and is an assumption that is made commonly for GRBs which turns out to be incorrect (see below). 

Thus, the first task should be the determination of the correlation between the variables. Unfortunately, 
most past works dealing with determination of the GRB LF or GRBFR (those mentioned above and others; 
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Fig. 1.— K-corrected Luminosity vs. redshift. The solid and dotted (black) curve shows the truncation due 
to a flux limits of /n m = 2 x 10~ 8 and 2 x 10 -9 erg/(s cm 2 ), respectively. In our analysis we use the larger 
and more conservative limit. The dashed (green) line shows the best fit luminosity evolution to the raw data 
(data points above the solid curve. Most of this correlation is due to the truncation. The solid (red) vertical 
and horizontal lines define the boundaries of the associated set (in this case M,) of the source marked by 
the red circle. The dotted (blue) lines and letters show the luminosity and redshift of the source and the 
maximum redshift and the minimum luminosity, defined in Eq. d6j) , that this source can have and still be in 
the sample. 


e.g. Kistler et al. 2007; Li 2008) omit this crucial step, which, as described below can lead to a incorrect 
GRBFR. A commonly used non-parametric method for determining correlations is the Spearman Rank test. 
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This method, however, fails for truncated data. Efron & Petrosian (1992) developed a novel method to 
account for this truncation. The gist of the method is to determine the rank i?j of each data point Li (or zi) 
among the Ni (or Mi) members of its associated set , or its largest un-truncated subset. The region bounded 
by the red lines in Fig. |T] depicts the boundary of this set for the data set L i: Zi identified by the red circle. 
This is for ranking in redshifts (from low to high values) and the set includes Mi sources with Zj < Zi and 
2 m axj > Zi (or Lj > Emin, i). In an analogous manner one can define the set of Ni sources with Lj > Li and 
L m 'i-n,j{z) < Li (or Zj < z max ,j for ranking the luminosities from above. These limiting values, also shown in 
Fig. m are given as 

-bmin J = 4/rd j ( Z-j ) f\ i Tn / A (.Zj) and Lj (z) = 4 ’zdj. (z max ,j ) film / N (Zmax.j ) ■ (6) 

Then using a test statistic, e.g. Kendell’s tau defined as r = Yii(Rj — Ej)/\jT^Vj, one can determine the 
degree of correlation. Here Ej and V 3 are the expectation and variance of the ranks. A small value (r -C 1) 
would imply independence, and r > 1 would imply significant correlation. In the latter case one can redefine 
the variables, e.g. define a “local” luminosity L' = L/g{z) [with <?(0) = 1] using a parametric form, e.g. 
g(z) = (1 + z) k , and calculate r as a function of the parameter fc . 4 The values of k for which r = 0 and 
r = ±1 give the best value and one sigma range for independence. 


Once independence is established then one can use the above mentioned non-parametric methods 
[(V/V max) or the C ~) to determine the mono-variate distributions; local LF ip(L') and the formation rate 
evolution p(z). We will use the C~ method which builds up the cumulative distributions of both variables 
L and z defined as 

(r{z) = J —^r dV d [ Z , ^ dz' and c/)(L') = J ii{x)dx (7) 

point by point non-parametrically, again using the concept of the associated set. For example 4>{Li) = 
n;= 2 (i + l/Nj). The cumulative functions can then be differentiated to get t/> and p using the derivative of 
the V(z), the co-moving volume up to z. 


The observed redshift distribution is related to the cumulative functions as 


dN 

dz 


dd{z) 

dz 


^minfc)] with L' min (z) = 


4?r d 2 L (z)f mil 
K{z)g{z) 


( 8 ) 


where (j)[L' min (z)\ represents the effects of the truncation and dcr(z)/dz gives the true redshift distribution of 
the parent population; i.e. what one would observe in the absence of truncation (i.e. when flux thresholds 

/lim 0). 


It should also be noted that, this method is not limited to a simple flux limited data but, as shown 
in Efron & Petrosian (1999), it can deal with truncations from both below and above, or with the most 
general truncation situation where each data point, say [Li,zi\, has its individual upper and lower limits, 
[Li t min < Li < Li t max and Zj, m in < Zj < Zj, max ]. We will not be dealing with such complications here, since 
our data is truncated from below only. 5 


In appendix A we apply this method to a flux limited sample selected from a parent simulated sample 
with known charcteristics of the LF and evolution, and demonstrate that we recover the input characteristics 
accurately. 


4 Note that for each k the limiting values L m jnj(z m axj) and associated sets Nj(Mj) are calculated anew. 

5 We note, however, that this is an enormous advantage because it allows one to combine data from many different regions 
of sky with different backgrounds and even from different instruments obtained with different selection criteria. 
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Fig. 2.— Variation of Kendell’s r statistics with the evolution index <5 showing independence for S = 2.33 
with one sigma range 1.48 to 3.12. We set S' = S and Z CI = 3.5 in Eq. (|9|). 


4. Results 


4.1. Luminosity Evolution 


As stressed above, the first task is to carry out test of the independence of L and z and determine the 
luminosity evolution. As shown above the raw data shows strong correlation L oc Z 3 9 . Using the associated 
sets to account for the truncations we find a Kendell’s r ~ 3.5 indicating the presence of a strong intrinsic 
correlation (i.e. luminosity evolution). The form of the evolution can be determined parametrically. The 
form often used is g(z ) = Z s which is appropriate for low and intermediate redshifts. But at higher redshifts 
this becomes excessively large considering the fact that the time intervals dt = dZ/(Z^/^I mZ 3 + 1 — SIm) 
decreases rapidly with increasing redshifts for Z > 3. Based on our experience in treating the evolution of 
the quasars we use the following, that has a slower luminosity evolution for Z > Z CI 


9(z) 


Z S x Zj 

Z* + ' 


( 9 ) 


After some trial and error we settled on Z cv = 3.5 and S' = 5, leaving us with one free parameter. We then 
define the “local” luminosity L' = L/g{z) (and the truncation curve L' min (z ) = L m i n (z)/g(z)) and calculate 
the statistic r as a function of S. Fig. [2] shows the results indicating independence of L' and z for S = 2.3 
with the one sigma range [1.5 — 3.1]. This is significantly less than the raw index ~ 4, but still considerable 
evolution; factors of 4.0, 7.8 and ~ 11 at redshifts 1, 2 and 3, respectively. This is in agreement with recent 











results using FF methods (e.g. Salvaterra et al. 2012) and earlier results using the Efron-Petrosian method 
on pseudo-redshift samples (Lloyd-Ronning et al. 2002; Yonetoku et al. 2004; Kocevski & Liang 2006) and 
samples with redshift (Petrosian et al. 2013). 


4.2. The Luminosity Function 


Transferring the observed L — z scatter diagram to the local luminosity L' — z diagram we can now 
determine the local LF ip{L') using the C~ method. Since we have assumed that the shape of the LF is 
invariant this function when shifted in luminosity by g(z ) describes the LF at all redshifts. The left panel 
of Fig. [3] shows the cumulative local LF 4>(L') build up point by point starting with the highest observed 
luminosity (filled points). For comparison we also show the raw cumulative distribution (open circles) 
showing the correction due to truncation obtained by our method. As shown in Appendix B one can also 
obtain a histogram of the differential LF 


#(£/) _ 4>(L') dlog <t>(L') 

dL' L' d log L' 


( 10 ) 


directly from the data. However, because of the paucity of the data and the random nature of the luminosities, 
their separation (e.g 5 log Li = log(Lj_i/Lj)) has a large dispersion yielding a noisy differential LF. Instead 
we obtain the differential LF by taking the derivative of appropriately smoothed curve fitted to </>(L'). As 
shown by the dashed (red) and dotted (green) curves in Fig. [3] (left) a smooth broken power law or a power 
law with exponential cutoff provide adequate fit to the cumulative LF. Derivatives of this give a similar 
description for the differential LF. For example the indicies -0.5 and -2.2 of </>(L') means that the differential 
luminosity function can also be fitted to a broken law with indicies dlogip/dlogL' ~ —1.5 and —3.2. This 
is a steeper LF, especially at the high end, than some recent results based on FF method; e.g. Howell et al. 
(2014) obtain indexes [—0.95,-2.59], and Salvattera et al. (2012) obtain indexes [—1.5,—2.3]. 


4.3. Formation Rate Evolution 


From the L' — z scatter diagram we can also obtain the GRBFR evolution. The cumulative rate 
evolution &(z) is shown by the filled points on the right panel of Fig. [3] This is obtained point by point 
starting with the lowest redshift. Again, for comparison we also show (open points) the raw observed 
cumulative distribution N(< z) = f~(dN/dz)dz. On the left panel of Fig. [I] we compare the smoothed 
differential observed redshift distribution dN/dz (dotted red) with the true distribution d&(z)/dz (solid, 
black), obtained from the derivative of a smooth solid (black) curve fitted to &{z) on Fig. [3] (right). Both 
of the above comparisons show that observations miss many intermediate and high redshift sources due 
to truncations affected by the gamma-ray threshold. These curves are normalized at low redshifts where 
the effect of truncation is expected to be minimal. A proper normalization can be carried out using total 
population counts like the so-called log-A-logS 1 diagram, which is beyond the scope of his paper. 


The co-moving GRBFR evolution is obtained as 


p{z) = Z 


da(z)/dz 
dV (. z)/dz 


( 11 ) 


This is shown by the solid (black) curve on the right panel of Fig. Q] In this figure we also show the raw 
GRBFR that one would obtain ignoring the truncation (dotted, red), i.e. using the observed dN/dz instead 
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Log Luminosity in erg/s 


No Luminosity Evolution__ 



Fig. 3.- Left: The cumulative LF. The filled points shows the cumulative local LF while the open 

points shows the cumulative distribution obtained ignoring the effects of the truncation. We have omitted 5 
GRBs with L < 10 49 erg/s. The dashed (red) curve shows a broken power law fit with indicies -0.5 and -2.2 
above and below Lq = 10 51 erg/s and the dotted (blue) curve shows a fit to Schechter function; a power-law 
with exponential cutoff. Right: The cumulative rate evolution a(z) versus redshifts corrected for selection 
effects due to the 7 -ray flux limit: filled points with solid (black) line a smoothed fit. The very similar 
dashed (green) curve is obtained using both gamma and X-ray threshold effects. The long-dashed (blue) is 
obtained ignoring luminosity evolution. The open points, and fitted dotted (red) curve, shows the observed 
cumulative distribution N(< z ), which of course does not include effects of the truncation. 


of da(z)/dz in Eq. (fill) . As expected this gives much lower rate at higher redshifts. We have also calculated 
GRBFR evolution ignoring the luminosity evolution, i.e. determining the distributions from the L — z data 
shown in Fig. [j] instead of L' — z scatter diagram. This is shown by the dashed (blue) curves in Figs. [3] 
and SI As evident ignoring the luminosity evolution overestimates the required GRBFR. This is the reason 
that some of earlier works (cited above) that ignored luminosity evolution obtained a high GRBFR at high 
redshifts, higher than the observed SFR. Again we have normalized these rates at low redshifts (z < 1) where 
the observational selection effects are smallest. 

All of the above results are obtained including only the effects truncation due to prompt 7 -ray detection 
threshold. As mentioned in § 2 , other observational selection effects, arising in the process of localization and 
securing redshifts, add more truncations. In particular the threshold of X-ray detection, the first crucial step 
for determination of the redshift, is an important factor. Based on the X-ray flux data from Racusin et al. 
(2011; private communications), we use a conservative limit of ~ 2 x 10 -13 erg/(s cm 2 ) (in the 0.2-10 keV 
range and at 11 hrs after the trigger) to determine how this additional selection bias affects our result. To 
this end we have repeated the above calculations using the effects of both gamma-ray and X-ray thresholds. 
This extra truncation redefines the associated sets which, as described above, for a single threshold consists 
of the set with Zj < Zi and z max j < Zi, where z max j, defined in Eq. J 6 |), is the maximum redshift the source 
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Fig. 4.— Left: Comparison of the observed redshift distribution dN/dz (dotted, red) with the intrinsic 
redshift distributions da/dz (solid, black) obtained from the differentiation of curves fitted to the cumula¬ 
tive distribution shown on the right panel of Fig. [3] The dashed (blue) curve shows the derived redshift 
distribution ignoring the luminosity evolution demonstrating the importance of inclusion of the luminosity 
evolution. The x points (green) are obtained including the effects of both gamma-ray and X-ray selection 
biases, as well of the luminosity evolution, as described in the text. This is very close to the black curve 
indicating that most of the selection bias is already accounted for by the correction for the gamma-ray 
threshold. These curves are normalized to redshifts z < 1 where the observational selection bias is expected 
to be small. Right: The density rate evolution obtained from the distributions on the left and Eq. ED using 
the same color codes. The dashed (red) curve is obtained using the smoothed observed dN/dz which shows 
a much more rapid decline than the true rate obtained for 7 -ray threshold only (solid black), both 7 and 
X-ray thresholds (long-dash green). The dashed (blue) curve is obtained ignoring the luminosity evolution. 
All rates are normalized at z < 1. 


with Lj , Zj can be moved to and still be in the sample. With two thresholds we now have two limiting 
redshifts. We use the minimum of the two to define the associated sets. The x points (green) on both panels 
of Figs. [I] show this results. These are very similar to the black curves obtained with the 7 -ray threshold 
effects alone. This is encouraging because it may imply that the other selection effects in securing a redshift, 
which cannot be easily quantified, have also small effects. Comparing these with the raw (uncorrected) 
redshift distributions (dotted red curve) we conclude that most of the correction has already been accounted 
for by the gamma-ray threshold. We also note that, there is a strong correlation between X-ray and optical 
fluxes of Swift GRBs which indicates that the X-ray threshold may be good proxy for the optical selection 
effects. 

Finally we note that a more rigorous approach to the multivariate nature of this problem would require 
us to obtain the evolution of the bivariate LF 4t(L 7 , Lx, z). As demonstrated in Singal et al. (2011) and 
(2013), the methods used in this paper for a single luminosity can be generalized to higher dimensions 
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without much difficulty. Similarly, if and when we understand the optical selection criterion, we can include 
optical luminosities as well in this multivariate distribution. 



Fig. 5.— Comparison of the the derived GRBFR (dash-dotted, green curve on Fig. [H right) and the SFR 
from Hopkins & Beacom (2006) normalized around redshift z ~ 3. The GRB rate (red) agrees remarkably 
well in the somewhat flat portion of the SFR rate and decays in a similar way at high redshifts. But at low 
redshifts it has a dramatically different shape and much higher rate. 


5. Summary and Discussion 

We have analyzed a sample of about 200 Swift flux limited long GRBs with known redshift using the 
non-parametric methods of Efron & Petrosian to obtain the general evolution of the luminosity function. 
Our results can be summarized as follows: 

• We find that the observed strong correlation between the luminosity and redshift is not totally due 
to the effects of the flux limit but that there is significant intrinsic correlation indicating that GRBs 
have undergone a strong luminosity evolution at least up to z < 3 which can be approximated as 
Lcx(l + z) 2 - 3 . 

• Correcting for this evolution we then obtain the local LF which can be represented by a broken power 
law or a Schechter type function. 

• We also find the true redshift distribution and the (co-moving) density rate evolution which starts at a 
maximum level at lowest redshifts, declines slowly up to 2 ~ 1 and merges into a plateau or an almost 
constant rate between redshift 1 and 3, and then shows a relatively steep decline at higher redshifts. 

• Thus, we find that GRBs undergo both luminosity and density evolution; they were more luminous but 
less numerous in the past than today. Often the luminosity evolution is ignored which as shown here 
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gives an incorrect (higher) density rate evolution at high redshifts. This is the source of many claims 
that GRBFR is higher than SFR at higher redshifts (see, e.g. Kistler et al. 2008; Salvaterra et al. 

2012; Wei et al. 2014). Robertson & Ellis (2012), and Dainotti et al. (2015) using afterglow emission, 
find a slower GRBR evolution, but their result also disagrees with ours because of their neglect of 
luminosity evolution. 

• Most analyses of GRB evolution assume a formation rate similar to the SFR especially at low redshift, 
but with an increase rate at higher z’s. The additional evolution is derived by FF to the observed 
distribution of redshifts and fluxes. With our non-parametric method we obtain directly the GRB 
formation rate which is quite different than those derived by FF methods. As shown in the right panel 
of Fig.0 our derived rate agrees very well with the SFR at redshifts z > 1, in particular in the plateau 
range 1 < 2 < 3, and shows a similar rapid decline at higher redshifts. But at low redshifts z < 1 
the GRB rate deviates by more than a decade from the SFR. This is puzzling because, according to 
the prevailing view that GRB formation is favored in low metalicity regions, one would expect a lower 
GRB rate in metal rich lower redshift galaxies. 

There are many potential factors that can account for this discrepancy. The above results are based on 
the LF of the prompt emission and includes rigorously only the y-ray threshold effects. But the data used 
involved also X-ray, UV and optical (space and ground based) measurement for localization and securing 
redshifts. We have included the X-ray detection threshold as best as can be done and find a small effect. The 
discrepancy at low redshifts remains large. It appears that the y-ray threshold effects account for most of the 
data truncations. This is encouraging because it may indicate that less quantifiable UV and optical selection 
biases may have similarly small effects. It should be noted that these added selection effects are expected to 
be more important at high redshifts so that they, most likely, are not the source of the discrepancy at low 
redshifts. The fact that all curves in Fig. [4] have same shape at low redshifts supports this view. 

Another possible factor may be the paucity of low redshift and low luminosity GRBs which dominate 
the shape of the nearby GRB rate. There have been some claims that low luminosity GRBs may belong to a 
separate class. However, as can be seen from Fig. [T] eliminating GRBs with L m j n < 10 49 erg/s from the list 
will affect only the result at z < 0.4 leaving a large discrepancy for 0.4 < z < 1. This discrepancy, therefore, 
may have an important implication on the formation of long GRBs. 

Acknowledgements: This work is partially supported by Swift guest investigator grants (NASA NNX12AE74G) 
and is based on Ellie Kitanidis’s senior honor thesis at Stanford University (see http: //purl. Stanford. edu/xp981bq5003). 


Appendix A: Testing the Methods 

To demonstrate the accuracy of our methods we have simulated a sample of sources with redshift and 
luminosity distributions similar to those of GRBs with specified form for the LF and density rate evolution 
but with no luminosity evolution. We have selected a flux limited subsample and applied the Efron-Petrosian 
methods to recover the intrinsic distributions. We recover no evolution (obtain Kendell’s r = 0 for k ~ 0.1), 
and as shown in Fig. [6l the calculated cumulative distributions 4>{L) and er(z), agree very well with the 
intrinsic distribution of the parent sample. 
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Fig. 6.— Comparison of the true intrinsic cumulative distributions (dashed)of a simulated sample (chosen 
to have characteristics similar to those of GRBs) with that of the flux limited “observed” sample (dotted) 
and that histogram obtained by our method that corrects for this selection bias. As evident the method 
reproduces the input forms extremely well, except at the extremities. Left: LF Right: density rate 

dota(z). 


Appendix B: Non-parametric Differential Distribution 


As mentioned in §3 the Efron-Petrosian method gives a non-parametric and point by point description 
of the cumulative distributions. Our usual procedure is to fit such histograms to a smooth function and take 
its derivative to obtain the differential distributions. This is because one cannot differentiate a discontinuous 
histogram. Here we describe a point by point estimation of the differential distributions (using luminosity 
as an example) which can give us the histogram of the differential distribution. 

The cumulative distribution up to and including each data point, say 4>{Li) is given as 


i=2 


( 12 ) 


where IV, is the number of data points in the associated set of point i (see Fig.[IJ). We can use two different 
estimator of the differential distribution il’(L) = —d(f){L)/dL. 

We can use the estimator = (4?(Lj) — 0(Lj-i)/(Li-i — L,) in which case it is easy to show that 

for Si = L i _ 1 /L l - 1 

<KLi) 1 


Liip(Li) = 

where the last equality is for JV* » 1. 


S t Ni (1 + 1 /Ni) SiNi 


-( 1-1 /Ni+ ...), 


(13) 


Another estimator is obtain using the fact that Eq. (3) yields <51n$(Lj) = In <f>{Li) — ln^»(Lj_i) = 
ln(l + 1/ATj), which then gives 

r in xt T M 1 + 1 / N i) _. <l>( L i) n ^ n a\ 

LMLi) - <t>{Li) dln ^ - <t>{Li) H1 + gi) - ^-(1 - 1/ (2JVi) - ^ + ..•), (14) 

which for A(> 1 and Si 1 is equal to the first estimator. 


We can similarly obtain the differential distribution da/dz and p(z). 


Unfortunately, for a sparse data set with sometimes large gaps, such as that available for GRBs, not 
all Ni are large (especially at the two ends of the data, and the relevant quantities In L — i or Si have large 
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dispersions. These fact introduce a large random noise so that we rely on smoothing and differentiation of 
the cumulative distributions. The above relation can be used for samples with dense and uniform coverage 
of the phase space. 
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